#include "OptData.h"

OptData::OptData(void) {
    this->m_nSizeX = 0;
    this->m_nSizeY = 0;
    this->m_nSizeZ = 0;
    this->m_nBeamletNumTotal = 0;
    //this->m_nWeightNumTotal=0;
}

OptData::OptData(const OptData &obj) {
    this->m_nSizeX = obj.m_nSizeX;
    this->m_nSizeY = obj.m_nSizeY;
    this->m_nSizeZ = obj.m_nSizeZ;
    this->m_nBeamletNumTotal = obj.m_nBeamletNumTotal;
    //this->m_nWeightNumTotal=obj.m_nWeightNumTotal;
    this->m_vFieldBeamletNum = obj.m_vFieldBeamletNum;
    this->m_vStructures = obj.m_vStructures;
    this->m_vDijIndex = obj.m_vDijIndex;
    this->m_vDijTotal = obj.m_vDijTotal;
    this->m_vWeightTotal = obj.m_vWeightTotal;
    this->m_vEUD = obj.m_vEUD;
}

OptData::~OptData(void) {
}

OptData& OptData::operator =(OptData& obj) {
    this->m_nSizeX = obj.m_nSizeX;
    this->m_nSizeY = obj.m_nSizeY;
    this->m_nSizeZ = obj.m_nSizeZ;
    this->m_nBeamletNumTotal = obj.m_nBeamletNumTotal;
    //this->m_nWeightNumTotal=obj.m_nWeightNumTotal;
    this->m_vFieldBeamletNum = obj.m_vFieldBeamletNum;
    this->m_vStructures = obj.m_vStructures;
    this->m_vDijIndex = obj.m_vDijIndex;
    this->m_vDijTotal = obj.m_vDijTotal;
    this->m_vWeightTotal = obj.m_vWeightTotal;
    this->m_vEUD = obj.m_vEUD;
    return *this;
}

void OptData::calculateDose(vector<double> *vWeight) {
    vector<vector<double> > *vDij = this->getDijTotal();
#ifdef USE_OMP
#pragma omp parallel for
#endif
    for (int i = 0; i < m_vDose.size(); i++) {
        m_vDose[i] = 0;
    }

    this->m_vDose.resize(vDij->size(), 0);

    int indexPoint = 0;
    int indexBeamlet;
#ifdef USE_OMP
#pragma omp parallel for default(shared) private(indexPoint,indexBeamlet) schedule(dynamic)
#endif
    for (indexPoint = 0; indexPoint < vDij->size(); indexPoint++) {
        double contriTemp = 0;
        double weiTemp = 0;
        for (indexBeamlet = 0; indexBeamlet < (*vDij)[indexPoint].size(); indexBeamlet++) {
            contriTemp = (*vDij)[indexPoint][indexBeamlet];

            weiTemp = (*vWeight)[indexBeamlet];

            //if contribution = 0 or weight =0 do not calculate
            if (!(((contriTemp >= -0.00000001) && (contriTemp <= 0.00000001)) || ((weiTemp >= -0.00000001) && (weiTemp <= 0.00000001)))) {
                this->m_vDose[indexPoint] += contriTemp*weiTemp;
            }
        }
    }
}

void OptData::calculateEUD(void) {
    if (m_vDose.size() == 0 || m_vStructures.size() == 0) {
        return;
    }
    for (int index = 0; index < m_vStructures.size(); index++) {
        double doseTotal = 0;
        int indexPoint = 0;
#ifdef USE_OMP
#pragma omp parallel for default(shared) reduction(+:doseTotal) private(indexPoint)
#endif
        for (indexPoint = 0; indexPoint < m_vStructures[index].getPoints()->size(); indexPoint++) {
            doseTotal += pow(m_vDose[(*m_vStructures[index].getPoints())[indexPoint]], m_vStructures[index].getEUDfac());
        }
        double eud = pow(doseTotal / (*m_vStructures[index].getPoints()).size(), 1 / m_vStructures[index].getEUDfac());
        this->m_vEUD.push_back(eud);
    }
}

vector<double> OptData::calculateDVH(vector<double> *resultDose, CStructure structure, double dvhResolution) {
    double maxDose = 0;
    for (int i = 0; i < resultDose->size(); i++) {
        if (maxDose < (*resultDose)[i]) {
            maxDose = (*resultDose)[i];
        }
    }
    int size = (int) (maxDose / dvhResolution + 1);

    vector<double> result;
    result.resize(size);

    for (int j = 0; j < structure.getPoints()->size(); j++) {
        double dose = (*resultDose)[(*structure.getPoints())[j]];
        int doseT = (int) (dose / dvhResolution);
        for (int index = 0; index <= doseT; index++) {
            result[index] = result[index] + 1;
        }
    }
    for (int j = 0; j < result.size(); j++) {

        result[j] = result[j] / structure.getPoints()->size();

    }
    return result;
}

void OptData::inputDoseSize(string inDoseSizePath) {
    pathCheck(inDoseSizePath);

    ifstream ifs;
    string line;

    try {

        ifs.open(inDoseSizePath.c_str());
        if (ifs.is_open()) {
            getline(ifs, line);
            if (line.length() > 0) {
                std::replace(line.begin(), line.end(), '.', ',');
                vector<string> vstrings = stringSplit(line, ",");
                if (vstrings.size() == 3) {
                    this->setDoseMatrixSize(
                            atoi(vstrings[0].c_str()), //x
                            atoi(vstrings[1].c_str()), //y
                            atoi(vstrings[2].c_str()) //z
                            );
                }
            }
        } else {
            throw runtime_error(inDoseSizePath + " file open error");
        }
        ifs.close();
        ifs.clear();
    } catch (exception ex) {
        ifs.close();
        ifs.clear();
        cerr << ex.what() << endl;
    }
}

void OptData::inputFieldType1(int fieldID, string inSpotFieldPath, string inSpotWeightPath) {

    pathCheck(inSpotFieldPath);
    pathCheck(inSpotWeightPath);

    ifstream ifsFieldFile;
    ifstream ifsWeightFile;
    string line;

    if ((this->m_vDijTotal.size()) != (m_nSizeX * m_nSizeY * m_nSizeZ)) {
        this->m_vDijTotal.resize(m_nSizeX * m_nSizeY * m_nSizeZ);
    }
    int baseBeamletNum = this->m_nBeamletNumTotal;
    int beamletNumField = 0;
    //read Weight File
    ifsWeightFile.open(inSpotWeightPath.c_str());
    if (ifsWeightFile.is_open()) {
        getline(ifsWeightFile, line);
        beamletNumField = atoi(line.c_str());
        this->m_vFieldBeamletNum.push_back(beamletNumField);
        while (getline(ifsWeightFile, line)) {
            if (line.length() > 0) {
                vector<string> vsWeightLine = stringSplit(line, ",");
                if (vsWeightLine.size() == 2) {
                    this->m_vWeightTotal.push_back(atof(vsWeightLine[1].c_str()));
                } else {
                    throw runtime_error("weight file error!");
                }
            }
        }
    } else {
        throw runtime_error(inSpotWeightPath + " file open error");
    }

    this->m_nBeamletNumTotal += beamletNumField;

    //for(vector< vector <double> >::size_type indexVDij = 0; indexVDij < m_vDijTotal.size(); indexVDij++)
    //{
    //	(m_vDijTotal)[indexVDij].resize(this->m_nBeamletNumTotal);
    //}

    try {
        unsigned int sizeX, sizeY, sizeZ;
        this->getDoseMatrixSize(sizeX, sizeY, sizeZ);
        if (sizeX == 0 || sizeY == 0 || sizeZ == 0) {
            throw runtime_error("Dose matrix size is needed to be read first!");
        }

        //read SpotField File

        ifsFieldFile.open(inSpotFieldPath.c_str());
        if (ifsFieldFile.is_open()) {
            getline(ifsFieldFile, line);
            while (getline(ifsFieldFile, line)) {
                if (line.length() > 0) {
                    vector<string> vsFieldLine = stringSplit(line, ",");
                    if (vsFieldLine.size() == 5) {
                        double contri = atof(vsFieldLine[4].c_str());

                        int point = this->position3DTo1D(atoi(vsFieldLine[0].c_str()),
                                atoi(vsFieldLine[1].c_str()), atoi(vsFieldLine[2].c_str()));

                        int beamletID = atoi(vsFieldLine[3].c_str());

                        if ((m_vDijTotal)[point].size() != this->m_nBeamletNumTotal) {
                            (m_vDijTotal)[point].resize(this->m_nBeamletNumTotal);
                        }
                        (m_vDijTotal)[point][beamletID + baseBeamletNum] = contri;
                    } else {
                        cout << line << endl;
                        throw runtime_error("field file error!");
                    }
                }
            }
        } else {
            throw runtime_error(inSpotFieldPath + " file open error");
        }

        ifsFieldFile.close();
        ifsWeightFile.close();

    } catch (exception ex) {
        ifsFieldFile.close();
        ifsFieldFile.clear();
        ifsWeightFile.close();
        ifsWeightFile.clear();

        cerr << ex.what() << endl;
    }
}

void OptData::inputFieldType2(int fieldID, string inSpotFieldPath, string inSpotWeightPath) {

    pathCheck(inSpotFieldPath);
    pathCheck(inSpotWeightPath);

    ifstream ifsFieldFile;
    ifstream ifsWeightFile;
    string line;

    if ((this->m_vDijTotal.size()) != (m_nSizeX * m_nSizeY * m_nSizeZ)) {
        this->m_vDijTotal.resize(m_nSizeX * m_nSizeY * m_nSizeZ);
    }
    int baseBeamletNum = this->m_nBeamletNumTotal;
    int beamletNumField = 0;
    //read Weight File
    ifsWeightFile.open(inSpotWeightPath.c_str());
    if (ifsWeightFile.is_open()) {
        getline(ifsWeightFile, line);
        beamletNumField = atoi(line.c_str());
        this->m_vFieldBeamletNum.push_back(beamletNumField);
        while (getline(ifsWeightFile, line)) {
            if (line.length() > 0) {
                vector<string> vsWeightLine = stringSplit(line, ",");
                if (vsWeightLine.size() == 2) {
                    this->m_vWeightTotal.push_back(atof(vsWeightLine[1].c_str()));
                } else {
                    throw runtime_error("weight file error!");
                }
            }
        }
    } else {
        throw runtime_error(inSpotWeightPath + " file open error");
    }
    //cout<<baseBeamletNum<<endl;
    this->m_nBeamletNumTotal += beamletNumField;
    //cout<<m_nBeamletNumTotal<<endl;

    //for(vector< vector <double> >::size_type indexVDij = 0; indexVDij < m_vDijTotal.size(); indexVDij++)
    //{
    //cout<<indexVDij<<" "<<this->m_nBeamletNumTotal<<endl;
    //	(m_vDijTotal)[indexVDij].resize(this->m_nBeamletNumTotal);
    //}

    try {
        unsigned int sizeX, sizeY, sizeZ;
        this->getDoseMatrixSize(sizeX, sizeY, sizeZ);
        if (sizeX == 0 || sizeY == 0 || sizeZ == 0) {
            throw runtime_error("Dose matrix size is needed to be read first!");
        }

        //read SpotField File
        ifsFieldFile.open(inSpotFieldPath.c_str());
        if (ifsFieldFile.is_open()) {
            while (getline(ifsFieldFile, line)) {
                if (line.length() > 0) {
                    vector<string> vsFieldLine = stringSplit(line, " ");
                    if (vsFieldLine.size() == 2) {
                        double contri = atof(vsFieldLine[1].c_str());

                        vector<string> vsFieldSubLine = stringSplit(vsFieldLine[0], ".");
                        if (vsFieldSubLine.size() == 5) {
                            int point = this->position3DTo1D(atoi(vsFieldSubLine[0].c_str()),
                                    atoi(vsFieldSubLine[1].c_str()), atoi(vsFieldSubLine[2].c_str()));

                            int beamletID = atoi(vsFieldSubLine[4].c_str());

                            if ((m_vDijTotal)[point].size() != this->m_nBeamletNumTotal) {
                                (m_vDijTotal)[point].resize(this->m_nBeamletNumTotal);
                            }
                            (m_vDijTotal)[point][beamletID + baseBeamletNum] = contri;
                        } else {
                            throw runtime_error("field file error!");
                        }
                    } else {
                        throw runtime_error("field file error!");
                    }
                }
            }
        } else {
            throw runtime_error(inSpotFieldPath + " file open error");
        }

        ifsFieldFile.close();
        ifsWeightFile.close();

    } catch (exception ex) {
        ifsFieldFile.close();
        ifsFieldFile.clear();
        ifsWeightFile.close();
        ifsWeightFile.clear();

        cerr << ex.what() << endl;
    }
}

bool OptData::inputFieldFromBin(int fieldID, string inSpotFieldPath, string inSpotWeightPath) {
    pathCheck(inSpotFieldPath);
    pathCheck(inSpotWeightPath);

    ifstream ifsFieldFile;
    ifstream ifsWeightFile;
    string line;

    if ((this->m_vDijTotal.size()) != (m_nSizeX * m_nSizeY * m_nSizeZ)) {
        this->m_vDijTotal.resize(m_nSizeX * m_nSizeY * m_nSizeZ);
    }
    int baseBeamletNum = this->m_nBeamletNumTotal;
    int beamletNumField = 0;
    //read Weight File
    try {
        ifsWeightFile.open(inSpotWeightPath.c_str());
        if (ifsWeightFile.is_open()) {
            getline(ifsWeightFile, line);
            beamletNumField = atoi(line.c_str());
            this->m_vFieldBeamletNum.push_back(beamletNumField);
            while (getline(ifsWeightFile, line)) {
                if (line.length() > 0) {
                    vector<string> vsWeightLine = stringSplit(line, ",");
                    if (vsWeightLine.size() == 2) {
                        this->m_vWeightTotal.push_back(atof(vsWeightLine[1].c_str()));
                    } else {
                        throw runtime_error("weight file error!");
                    }
                }
            }
        } else {
            throw runtime_error(inSpotWeightPath + " file open error");
        }

        this->m_nBeamletNumTotal += beamletNumField;

        //read Contri File
        unsigned int sizeX, sizeY, sizeZ;
        this->getDoseMatrixSize(sizeX, sizeY, sizeZ);
        if (sizeX == 0 || sizeY == 0 || sizeZ == 0) {
            throw runtime_error("Dose matrix size is needed to be read first!");
        }

        //read SpotField binary File
        ifsFieldFile.open(inSpotFieldPath.c_str(), ios::in | ios::binary);
        if (!ifsFieldFile.is_open()) {
            cerr << "ERROR opening the spot contribution file to read: " << inSpotFieldPath << endl;
        }

        int nx, ny, nz, nSpotCachedPerVox;
        ifsFieldFile.read((char*) &nx, sizeof (int));
        ifsFieldFile.read((char*) &ny, sizeof (int));
        ifsFieldFile.read((char*) &nz, sizeof (int));
        if (sizeX != nx || sizeY != ny || sizeZ != nz) {
            cerr << "Error: The binary cache file does not match the Scan Ion Control Point Sequence" << endl;
            return false;
        }

        int nRCursor = 0;
        while (!ifsFieldFile.eof() && nRCursor < sizeX * sizeY * sizeZ) {
            ifsFieldFile.read((char*) &nx, sizeof (int));
            ifsFieldFile.read((char*) &ny, sizeof (int));
            ifsFieldFile.read((char*) &nz, sizeof (int));
            if (nx >= sizeX || ny >= sizeY || nz >= sizeZ) {
                cerr << "Error: the dose voxel index input from the file is out of range: (" << nx << ", " << ny << ", " << nz << ")" << endl;
                //return false;
                continue;
            }

            ifsFieldFile.read((char*) &nSpotCachedPerVox, sizeof (int));
            if (nSpotCachedPerVox <= 0) continue;

            int point = this->position3DTo1D(nx, ny, nz);
            if ((m_vDijTotal)[point].size() != this->m_nBeamletNumTotal) {
                (m_vDijTotal)[point].resize(this->m_nBeamletNumTotal);
            }
            INFLU_ELEM *arrInfluElems = new INFLU_ELEM[nSpotCachedPerVox];
            ifsFieldFile.read(reinterpret_cast<char*> (arrInfluElems), sizeof (INFLU_ELEM) * nSpotCachedPerVox);

            for (int ispt = 0; ispt < nSpotCachedPerVox; ispt++) {
                long nSpotID = arrInfluElems[ispt].m_nID;
                float fContrib = arrInfluElems[ispt].m_fContrib;

                if ((nSpotID)>this->m_nBeamletNumTotal) {
                    cerr << "Error: the beamlet index input from the file is out of range: (" << nSpotID << ")" << " Total = " << this->m_nBeamletNumTotal << endl;
                    cerr << nRCursor << endl;
                    //return false;
                    continue;
                }

                (m_vDijTotal)[point][nSpotID] = fContrib;
            }

            delete [] arrInfluElems;
            nRCursor++;
        }
        ifsFieldFile.close();
        ifsWeightFile.close();

    } catch (exception ex) {
        ifsFieldFile.close();
        ifsFieldFile.clear();
        ifsWeightFile.close();
        ifsWeightFile.clear();

        cerr << ex.what() << endl;
    }
}

void OptData::inputField(int fieldID, string inSpotFieldPath, string inSpotWeightPath) {
    pathCheck(inSpotFieldPath);

    ifstream ifsFieldFile;
    string line;
    ifsFieldFile.open(inSpotFieldPath.c_str());
    if (int pos = inSpotFieldPath.find("cabin") != string::npos) {
        inputFieldFromBin(fieldID, inSpotFieldPath, inSpotWeightPath);
    } else if (ifsFieldFile.is_open()) {
        getline(ifsFieldFile, line);
        if (line.length() > 0) {
            vector<string> vsFieldLine = stringSplit(line, ",");
            if (vsFieldLine.size() == 3) {
                inputFieldType1(fieldID, inSpotFieldPath, inSpotWeightPath);
            } else {
                inputFieldType2(fieldID, inSpotFieldPath, inSpotWeightPath);
            }
        }
    } else {
        throw runtime_error(inSpotFieldPath + " file open error");
    }

    ifsFieldFile.close();
}

void OptData::inputStructure(string structName, string inStructPath, double EUD) {
    pathCheck(inStructPath);

    ifstream ifs;
    string line;

    CStructure pStruct;
    //vector<int> vPoints = pStruct->getPoints();
    pStruct.setStructName(structName);
    pStruct.setEUDfactor(EUD);

    unsigned int x, y, z;
    this->getDoseMatrixSize(x, y, z);
    try {

        ifs.open(inStructPath.c_str());
        if (ifs.is_open()) {
            while (getline(ifs, line)) {
                if (line.length() > 0) {
                    std::replace(line.begin(), line.end(), '.', ',');
                    vector<string> vstrings = stringSplit(line, ",");
                    if (vstrings.size() == 3) {

                        int positionTemp = 0;

                        if ((atoi(vstrings[0].c_str()) < x) && (atoi(vstrings[1].c_str()) < y) && (atoi(vstrings[2].c_str()) < z)) {
                            positionTemp = this->position3DTo1D(atoi(vstrings[0].c_str()),
                                    atoi(vstrings[1].c_str()), atoi(vstrings[2].c_str()));

                            pStruct.getPoints()->push_back(positionTemp);
                            //vPoints->push_back(positionTemp);
                        }
                    }
                }
            }
        } else {
            throw runtime_error(inStructPath + " file open error");
        }

        // add structure to matrix
        this->addStructure(pStruct);

        //delete pStruct;
        ifs.close();
        ifs.clear();
    } catch (exception ex) {
        ifs.close();
        ifs.clear();

        //delete pStruct;
        cerr << ex.what() << endl;
    }
}

void OptData::outputWeight(vector<double> *vWeight, string outputDir) {

    path pathOutput(outputDir);
    try {
        if (!exists(pathOutput)) {
            create_directory(pathOutput);
        }
        if (exists(pathOutput)) // does p actually exist?
        {
            if (is_directory(pathOutput)) // is p a directory?
            {
                ofstream ofs;
                for (int indexOfField = 0; indexOfField<this->m_vFieldBeamletNum.size(); indexOfField++) {
                    //calculate offset
                    int offset = 0;
                    for (int indexOffset = 0; indexOffset < indexOfField; indexOffset++) {
                        offset += m_vFieldBeamletNum[indexOffset];
                    }
                    int endPoint = m_vFieldBeamletNum[indexOfField] + offset;

                    std::stringstream pathBuff;
                    //pathBuff.str("");
                    pathBuff << "resultWeight_" << indexOfField + 1 << ".txt";
                    string strWtFullPath = (pathOutput / path(pathBuff.str())).string();
                    ofs.open(strWtFullPath.c_str());
                    //cout<<strWtFullPath.c_str()<<endl;
                    if (ofs.is_open()) {
                        ofs << endPoint - offset << endl;
                        for (int indexInTotalWeight = offset; indexInTotalWeight < endPoint; indexInTotalWeight++) {
                            ofs << indexInTotalWeight - offset << "," << (*vWeight)[indexInTotalWeight] << endl;
                        }
                    }

                    ofs.close();
                }
            } else {
                cout << pathOutput << " exists, but is neither a regular file nor a directory\n";
            }
        } else
            cout << pathOutput << " does not exist\n";
    } catch (const filesystem_error& ex) {
        cout << ex.what() << endl;
    }
}

void OptData::outputDij(vector<vector<double> > *dij, string outDijPath) {


    ofstream ofs;

    ofs.open(outDijPath.c_str());
    if (ofs.is_open()) {
        for (int indexVDij = 0; indexVDij != dij->size(); indexVDij++) {
            for (int indexDij = 0; indexDij != (*dij)[indexVDij].size(); indexDij++) {
                ofs << (*dij)[indexVDij][indexDij] << " ";
            }
            ofs << endl;
        }
    }

    ofs.close();
    ofs.clear();
}

void OptData::outputDose(vector<double> *dose, string outputDir) {

    path pathOutput(outputDir);
    try {
        if (!exists(pathOutput)) {
            create_directory(pathOutput);
        }
        if (exists(pathOutput)) // does p actually exist?
        {
            if (is_directory(pathOutput)) // is p a directory?
            {
                ofstream ofs;

                std::stringstream pathBuff;
                //pathBuff.str("");
                pathBuff << "resultDose.txt";
                string strDoseFullPath = (pathOutput / path(pathBuff.str())).string();
                ofs.open(strDoseFullPath.c_str());
                //cout<<strWtFullPath.c_str()<<endl;
                if (ofs.is_open()) {
                    ofs << endl;
                    ofs << "#0" << endl;
                    for (int i = 0; i != dose->size(); i++) {
                        ofs << (*dose)[i] << " ";
                        if ((i != 0) && (i != dose->size() - 1) && ((i + 1) % m_nSizeX) == 0) {
                            ofs << endl;
                        }
                        if ((i != 0) && (i != dose->size() - 1) && ((i + 1) % (m_nSizeY * m_nSizeX)) == 0) {
                            ofs << endl;
                            ofs << "#" << ((i + 1) / (m_nSizeY * m_nSizeX)) << endl;
                        }
                    }
                }

                ofs.close();
                ofs.clear();
            } else {
                cout << pathOutput << " exists, but is neither a regular file nor a directory\n";
            }
        }
    } catch (const filesystem_error& ex) {
        cout << ex.what() << endl;
    }
}

void OptData::outputObjValue(vector<double> *objValue, string outputDir) {

    path pathOutput(outputDir);
    try {
        if (!exists(pathOutput)) {
            create_directory(pathOutput);
        }
        if (exists(pathOutput)) // does p actually exist?
        {
            if (is_directory(pathOutput)) // is p a directory?
            {
                ofstream ofs;

                std::stringstream pathBuff;
                //pathBuff.str("");
                pathBuff << "objValue.txt";
                string strDoseFullPath = (pathOutput / path(pathBuff.str())).string();
                ofs.open(strDoseFullPath.c_str());
                //cout<<strWtFullPath.c_str()<<endl;
                if (ofs.is_open()) {

                    for (int i = 0; i != objValue->size(); i++) {
                        ofs << (*objValue)[i] << endl;
                    }
                }

                ofs.close();
                ofs.clear();
            } else {
                cout << pathOutput << " exists, but is neither a regular file nor a directory\n";
            }
        }
    } catch (const filesystem_error& ex) {
        cout << ex.what() << endl;
    }
}

void OptData::outputDVH(vector<double> *resultDose, string outputDir, double dvhResolution) {
    path pathOutput(outputDir);
    try {
        if (!exists(pathOutput)) {
            create_directory(pathOutput);
        }

        if (exists(pathOutput)) // does p actually exist?
        {
            if (is_directory(pathOutput)) // is p a directory?
            {
                ofstream ofs;

                vector<vector<double> > dvhs;
                int maxDVHSize = 0;
                //pathBuff.str("");
                for (int index = 0; index < m_vStructures.size(); index++) {
                    std::stringstream pathBuff;
                    pathBuff << m_vStructures[index].getStructName() << ".dvh";
                    string strDoseFullPath = (pathOutput / path(pathBuff.str())).string();
                    ofs.open(strDoseFullPath.c_str());

                    vector<double> dvh = calculateDVH(resultDose, m_vStructures[index], dvhResolution);
                    cout << m_vStructures[index].getStructName() << " " << dvh.size() << endl;

                    dvhs.push_back(dvh);
                    maxDVHSize = (maxDVHSize < dvh.size()) ? dvh.size() : maxDVHSize;


                    for (int i = 0; i < dvh.size(); i++) {
                        if (dvh[i] != 0) {
                            ofs << (double) i * dvhResolution << "\t" << dvh[i] << endl;
                        }
                    }
                    ofs.close();
                    ofs.clear();
                }

                cout << maxDVHSize << endl;
                for (int index = 0; index < m_vStructures.size(); index++) {
                    dvhs[index].resize(maxDVHSize, 0);
                }

                std::stringstream pathBuff;
                pathBuff << "total.dvh";
                string strDoseFullPath = (pathOutput / path(pathBuff.str())).string();
                ofs.open(strDoseFullPath.c_str());

                ofs << "dose" << "\t";
                for (int index = 0; index < m_vStructures.size(); index++) {
                    ofs << m_vStructures[index].getStructName() << "\t";
                }
                ofs << endl;

                for (int i = 0; i < maxDVHSize; i++) {
                    ofs << (double) i * dvhResolution << "\t";
                    for (int idvh = 0; idvh < dvhs.size(); idvh++) {
                        ofs << dvhs[idvh][i] << "\t";
                    }
                    ofs << endl;
                }
                ofs.close();
                ofs.clear();

            } else {
                cout << pathOutput << " exists, but is neither a regular file nor a directory\n";
            }
        }
    } catch (const filesystem_error& ex) {
        cout << ex.what() << endl;
    }
}

void OptData::inputDose(string inDosePath) {
    pathCheck(inDosePath);

    ifstream ifs;
    string line;

    try {

        ifs.open(inDosePath.c_str());
        if (ifs.is_open()) {
            while (getline(ifs, line)) {
                if (line.length() > 0) {
                    vector<string> vstrings = stringSplit(line, " ");
                    if (vstrings.size() > 1) {
                        for (int i = 0; i < vstrings.size(); i++) {
                            m_vDose.push_back(atof(vstrings[i].c_str()));
                        }
                    }

                }
            }
        } else {
            throw runtime_error(inDosePath + " file open error");
        }

        ifs.close();
        ifs.clear();
    } catch (exception ex) {
        ifs.close();
        ifs.clear();

        //delete pStruct;
        cerr << ex.what() << endl;
    }
}

